Discovery of inhibitors of protein tyrosine phosphatase 1B contained in a natural products library from Mexican medicinal plants and fungi using a combination of enzymatic and in silico methods**

This work aimed to discover protein tyrosine phosphatase 1B (PTP1B) inhibitors from a small molecule library of natural products (NPs) derived from selected Mexican medicinal plants and fungi to find new hits for developing antidiabetic drugs. The products showing similar IC50 values to ursolic acid (UA) (positive control, IC50 = 26.5) were considered hits. These compounds were canophyllol (1), 5-O-(β-D-glucopyranosyl)-7-methoxy-3′,4′-dihydroxy-4-phenylcoumarin (2), 3,4-dimethoxy-2,5-phenanthrenediol (3), masticadienonic acid (4), 4′,5,6-trihydroxy-3′,7-dimethoxyflavone (5), E/Z vermelhotin (6), tajixanthone hydrate (7), quercetin-3-O-(6″-benzoyl)-β-D-galactoside (8), lichexanthone (9), melianodiol (10), and confusarin (11). According to the double-reciprocal plots, 1 was a non-competitive inhibitor, 3 a mixed-type, and 6 competitive. The chemical space analysis of the hits (IC50 < 100 μM) and compounds possessing activity (IC50 in the range of 100–1,000 μM) with the BIOFACQUIM library indicated that the active molecules are chemically diverse, covering most of the known Mexican NPs’ chemical space. Finally, a structure–activity similarity (SAS) map was built using the Tanimoto similarity index and PTP1B absolute inhibitory activity, which allows the identification of seven scaffold hops, namely, compounds 3, 5, 6, 7, 8, 9, and 11. Canophyllol (1), on the other hand, is a true analog of UA since it is an SAR continuous zone of the SAS map.


Introduction
More than 474 million people worldwide live with type 2 diabetes mellitus (T2DM), characterized by chronic hyperglycemia due to insulin resistance and pancreatic β-cell dysfunction.Furthermore, in 2021, the International Diabetes Federation has estimated that over 6 million deaths yearly are due to T2DM (IDF, 2021).To control high blood glucose levels in diabetic patients, in addition to a healthy lifestyle, drug treatments and/or insulin are prescribed.Nonetheless, in some countries such as in Mexico, a large segment of the population uses medicinal plants for treating diabetes, alone or in combination with allopathic medicines.
The global prevalence of diabetes has prompted the search for novel and more efficacious therapeutic agents.In this regard, a better understanding of the complicated mechanisms involved in T2DM has inspired the pursuit of new compounds targeting new receptors and crucial enzymes involved in glucose homeostasis (ADA, 2022).Protein tyrosine phosphatase 1B (PTP1B) shows promise among these enzymes.PTP1B is a critical enzyme in the dephosphorylation of the insulin receptor and its downstream signaling components.The relevance of this enzyme in the pathophysiology of diabetes has been demonstrated by the genetic deletion of PTP1B in mice; PTP1Bdeficient mice remain insulin sensitive on a high-fat diet when compared with the wild-type (Na et al., 2016;Kanwal et al., 2022).
PTP1B is composed of three domains, an N-terminal catalytic domain (1-300) that contains the crucial Cys215, Asp181, and Gln262 residues; a regulatory domain (301-400) that is prolinerich and confers substrate specificity to PTP1B; and finally, a C-terminal domain (401-435) that is involved in the binding to the plasmatic reticulum and is intrinsically disordered (Tonks, 2003;Song et al., 2017;Liu et al., 2022;Quy et al., 2022).Since the nonreceptor protein tyrosine phosphatase family has a highly conserved catalytic site (~75% of sequence similarity), searching for selective allosteric inhibitors is crucial in drug discovery.
To discover new hit compounds that inhibit PTP1B, screening libraries based on natural products, mostly from antidiabetic plants, is an excellent and straightforward strategy (Singh et al., 2022).In addition, natural product libraries occupy a greater area of chemical space than synthetic compounds, reduce screening costs, and increase the speed of finding appropriate hits (Lahlou, 2007;Harvey et al., 2010).Finally, it is important to point out that a few natural products and some chemical derivatives, such as trodusquemine from the dogfish shark, ertiprotafib, and JTT-551, have entered clinical trials (Kazakova et al., 2022).
Based on the aforementioned considerations, this work aimed to 1) identify molecules from a small molecule library of natural products obtained from Mexican medicinal plants and fungi that might inhibit PTP1B by in-house enzymatic screening; 2) explore their chemical space; and 3) identify the scaffold hops and activity with the structure-activity similarity (SAS) map.Altogether, all these ensure a high success rate of finding new PTP1B inhibitors with chemically original or similar structures and increase the probability of identifying hit compounds.

PTP1B inhibition assay
A spectrophotometric and colorimetric method was used to detect the inhibitory activity of 99 in-house compounds on hPTP1B 1-400 (Rangel-Grimaldo et al., 2020).The analyses were carried out in triplicate; the test materials and positive control (ursolic acid) were dissolved in DMSO and buffer solution.Ursolic acid was used as a reference since it has been widely reported as one of the best allosteric PTP1B inhibitor by different research groups (Liu et al., 2022;Quy et al., 2022).The compounds were initially tested in concentrations between 20 and 1,000 μM, enzyme (66 nM) and buffer [Tris-HCl 50 mM, pH 6.8] solution, and 0.125 mM of 4-nitrophenyl phosphate (4-NP, Sigma-Aldrich, St. Louis, MO, United States) was incubated at room temperature for 20 min as previously described.At the end of the incubation, the reaction was terminated by adding 5 μL of NaOH (10 M); then, the absorbance was measured at 405 nm (t20).The inhibitory activity was determined as a percentage compared to the blank (Buffer) according to the following equation: where A 405(s) is the ΔA of the sample (A t20 − A t0 ), A 405(b) is the ΔA of the blank obtained in the same manner as that of the sample.Using the data obtained, the enzyme inhibition curve was plotted, for which the average percentage of inhibition was plotted as a function of the inhibitor concentration.To calculate the IC 50 , non-linear regression was performed using the Origin 8.0 software (OriginLab Corporation, 2022).For calculating the IC 50 , the compounds showing a higher percent of inhibition at 20 μM were tested at different concentrations, as shown in Supplementary Figure S1.

Compound pIC 50
[a] ΔA [b]  Compound pIC 50 [a] ΔA [b]    Negative log of the medium inhibitory concentration (IC 50 ) in Molar. [b] Activity difference between UA and the compound.

In-house library of natural products
The tested in-house compounds (1-99) were isolated from medicinal plants and fungi, as mentioned in Supplementary Table S1, which includes the natural sources and pertinent references.These compounds belong to a small molecule library of natural products from the Department of Pharmacy, School of Chemistry, UNAM.

Preparation of natural products databases
The isomeric SMILES codes of the most active tested compounds were generated using the Open Babel toolbox (O'Boyle et al., 2011).For chemical space characterization, we included our in-house library (most active NPs) and 423 Mexican natural products from the BIOFACQUIM database (Pilón-Jiménez et al., 2019).

Molecular fingerprints
Six molecular fingerprints (FPs) were calculated for the natural products database employing the MayaChemTools (Sud, 2016), rcdk (Guha, 2007), and RDKit (Landrum, 2022) packages (Supplementary Figure S2).These FPs were further used to compute the Tanimoto similarity index from the similarity matrix analysis.The Tanimoto similarity cumulative distribution (TSD) plots were generated to select the best FP representation for the tested compounds and the BIOFACQUIM library (Supplementary Figure S2).The PubChem bit-string (881 bits) was computed for the natural products database chemical space representation using the Rcpi (Cao et al., 2015) module of the R software (R Core Team, 2020).
In the case of the most active NPs, 25 molecular FPs were calculated, and the same methodology was followed to select the best FP representation and plot the SAS map (Supplementary Figure S14).

Kinetic analysis
The conditions for preparing reagents (enzyme and substrate), incubation, and analysis of results were the same as that described in Section 2.1 (Jiménez-Arreola et al., 2020;Díaz-Rojas et al., 2021).Additional steps in the methodology are described below.
To characterize the type of inhibition exerted by 1, 3, and 6 on PTP1B, the hydrolysis of 4-NP at the sub-saturating state and increasing concentrations of each compound was measured.An enzyme saturation curve was prepared with a 4-NP stock solution (10 mM); then, a series of curves with variable 4-NP concentrations and at least five fixed concentrations of the inhibitors were built; by considering the IC 50 value as the midpoint, the type of inhibition of 1, 3, and 6 was estimated; each point of the curves was obtained in triplicate.The compound concentration yielding 50% activity inhibition was interpolated from the results by fitting into a general inhibition model.The resulting set of activities was fitted globally using the non-linear regression algorithm of Levenberg-Marquardt as implemented in the Gnuplot (Williams et al., 2017).Each group was fitted to the kinetic equations for linear competitive, linear uncompetitive, linear non-competitive (classic), linear mixed type, parabolic competitive, parabolic uncompetitive, and parabolic noncompetitive inhibition mechanisms.The best fit was ascertained by the reduced chi-squared (χ 2 ) value, residual distribution, and uncertainty in the parameter estimates.The best description of the inhibition data for compounds 1, 3, and 6 was obtained with Eq. 2, Eq. 3, and Eq. 4, respectively, where V MAX , K M , and K Icu , K IC , and K IU are the maximum velocity, Michaelis constant for the substrate 4-NP, and inhibition constants for each type of inhibition, respectively.

Structure-activity similarity
The SAS map of the most active compounds tested against PTP1B was constructed using the absolute value of the pIC 50  difference (|ΔActivity|) and the UA similarity calculated with the atom type fingerprints (Supplementary Figure S14), following our previous work methodology (Santiago et al., 2021).The SAS map can only be constructed using the absolute value of the pIC 50 difference (|Δ Activity|) against ursolic acid and the UA similarity of Tanimoto (Bajusz et al., 2015).

Structural model of PTP1B 1-400
A structural model of the hPTP1B 1-400 protein was obtained from the AlphaFold Protein Structure Database developed by DeepMind and EMBL-EBI (https://alphafold.ebi.ac.uk/).The UniProt code P18031 corresponds to the PTPN1 gene, which codes for human PTP1B that is among the proteins in the database predicted by AlphaFold (code: Q9PT91).The pdb file was downloaded from the following link: https://alphafold.ebi.ac.uk/entry/P18031 (David et al., 2020).The coordinates of this model were prepared to perform a molecular dynamics (MD) simulation using the LEAP module from AmberTools 2021.The structure was submitted to the following procedure: hydrogens were added using the LEAP module with the leaprc.protein.ff19SBforce field; K + counter ions were also included to neutralize the system.The protein was solvated in an octahedral box of explicit TIP3P model water molecules localizing the box limits at 12 Å from the protein surface.MD simulations were performed at 1 atm and 315 K, and maintained with the Berendsen barostat and thermostat, respectively, using periodic boundary conditions and particle mesh Ewald sums (grid spacing of 1 Å) for treating long-range electrostatic interactions with a 10-Å cutoff for computing direct interactions.The SHAKE algorithm was used to satisfy bond constraints, using a time step of two femtoseconds (fs) to integrate Newton's equations as recommended in the Amber package.All calculations were made using the graphics processing unit (GPU)accelerated MD engine in Amber (pmemd.cuda), a program package that runs entirely on CUDA-enabled GPUs (Case et al., 2005(Case et al., , 2012)).The protocol minimized the initial structure, followed by 50picosecond (ps) heating and pressure equilibration at 315 K and 1. 0 atm pressure, respectively.Finally, the system was equilibrated with 500 ps before starting the production of MD.The production of the MD consisted of 100 nanoseconds (ns).The frames were saved at 10ps intervals for subsequent analysis.All analyses were done using CPPTRAJ (Roe and Cheatham, 2013).Root mean square deviations (RMSDs) and root mean square fluctuations (RMSFs) were calculated, considering the C, Cα, and N. The charts were built using OriginPro 9.1, and the trends were adjusted with smooth function processing (lowess span method).

Compound
Reduce

Molecular docking
The docking analysis was done using the structural model of hPTP1B 1-400 obtained after MD simulation.Structures 1-11 were constructed and minimized using the Avogadro software (Hanwell et al., 2012).The AutoDockTools 1.5.4 was used to prepare the pdb files of the protein and compounds.Polar hydrogen atoms and Kollman united-atom partial charges were added to the protein structures.By contrast, Gasteiger-Marsili charges and rotatable groups were automatically assigned to the structures of the ligands.We used AutoDock Vina to carry out the docking, covering the entire enzyme.The grid box size was 126 Å × 126 Å × 126 Å in the x, y, and z dimensions and the central coordinates of 54.94, 68.66, and 63.45 for x, y, and z, respectively, with exhaustiveness of 8.The best conformational states were visualized by using PyMOL version 2.4.0 and Maestro version 5. 3.156 (DeLano, 2004).
In the case of compound 6, the geometric isomers and their more stable conformers (6a1, 6a2, 6b1, and 6b2) were optimized using SQM-MP7 methodology with implicit aqueous solvation and finally converted to PDBQT format.The three-dimensional model of PTP1B was the same as that used in the aforementioned compounds.A wide box comprised the whole putative active site and its neighborhood within 7 Å.Several docking rounds were performed using AutoDock Vina 1.2. 3 (Eberhardt et al., 2021) to collect at least 160 docking poses for each compound.The compounds were clustered using the Clustering Tool plugin for VMD (Humphrey et al., 1996), using a 1-Å cutoff, and the average Vina docking energy of each cluster was obtained from Vina logs.Each cluster's RMSD was calculated using one arbitrary pose as a reference.Clusters having one or more poses with docking energy higher than −5.3 kcal/mol were ignored because the more stable ligand-protein complexes for compound 6 presented more negative values, therefore all values above −5.3kcal/mol were disregarded.The docking events with lower energy fell into two neighboring but well-separated binding pockets designated as I and II.For each binding pocket, the conformation with lower docking energy and lower RMSD was considered more likely to represent the real conformation, although docking poses occupying each pocket were considered.

FIGURE 6
Chemical structure of the 19 pairs of compounds identified as activity cliffs in the SAS map.The compared pairs are represented by a gray line, and the molecular similarity value is highlighted in yellow.TABLE 3 Results of molecular docking between PTP1B and compounds 1, 3, and 6.

Acute oral toxicity in mice
The potential acute toxicity of compound 6 was determined according to the Lorke's method (Lorke, 1983).The evaluation was performed following the Official Mexican Norm for Laboratory Animal Care and Use (NOM-062-ZOO-1999), with international conventional codes for laboratory animal use, and approved by the Institutional Committee for Care and Use of Laboratory Animals (CICUAL-FQ), Facultad de Química, UNAM (FQ/CICUAL/440/ 21).Eight-week-old male ICR mice (25-36 g) were adjusted to laboratory conditions.At the end of the experiments, the mice were euthanized by hypoxia in a CO 2 chamber.
Compound 6 was fed by an intragastric route in two independent phases.In each case, 12 mice were divided into four groups (n = 3).In the first stage, the animals were treated with 10, 100, and 1,000 mg/kg; in the second phase, 1,600, 2,900, and 5,000 mg/kg were administered.The control animals were fed 0.05% Tween 80 ® in saline solution.The weight of the animals was measured daily for 14 days at each stage.At the end of the test, all animals were euthanized to obtain the lungs, heart, kidneys, and liver to detect macroscopic organ injury.The mice were observed to identify acute toxic effects, changes in behavioral patterns, or mortality.

Inhibition of PTP1B and chemical space of NPs
A screening using hPTP1B 1-400 of a small customized (inhouse) NP library containing 99 compounds, mostly from medicinal plants, was performed.The molecules were selected by the following criteria: ethnopharmacological background of the plant source, good physicochemical properties such as solubility, and the availability of the compounds for future research.The preliminary screening process at two different concentrations (20 and 1,000 μM) allowed the identification of those molecules that showed 90% inhibition or more at 20 μM.Thereafter, the IC 50 values of these molecules were calculated.The results allowed us to identify 47 inhibitory NPs, with IC 50 values ranging from 30 to 1,000 μM, and the most active were compounds 1-11, whose activity gap ranged between 30 and 100 μM.Ursolic acid was used as the reference since it has been widely reported as one of the best allosteric PTP1B inhibitor by different research groups (Liu et al., 2022;Quy et al., 2022).Next, their pIC 50 values [the negative log of the medium inhibitory concentration (IC 50 )] (Table 1, Supplementary Table S2, and Supplementary Figure S1) were calculated; the pIC 50 value is a more reliable parameter to compare the potency of the different compounds tested at the same molar levels, and it is commonly used in chemoinformatics studies; furthermore, the pIC 50 value covers both micro-and millimolar potencies (Kalliokoski et al., 2013;Thakur et al., 2022).

Kinetic analysis of selected compounds
Compounds 1, 3, and 6 were selected for kinetic analysis due to their sufficient availability and good inhibitory activity against PTP1B (Table 1; Table 2; Figure 3; Figure 4).Compound 1 (the most active compound identified, according to pIC 50 and ΔA) behaved as a non-competitive inhibitor and 3 as a mixed inhibitor (with the second-best value of pIC 50 and ΔA).In both cases, when the concentration of the compounds increased, the slope and intersection with the ordinate axis changed, but the intersection with the abscissae did not.In addition, the K IC and K IU values changed for compound 3 (Spector and Cleland, 1981;Copeland, 2000).In both cases, the inhibition mechanism discovered was relevant because they interacted with a distinct area from the active site, guaranteeing a selective interaction with PTP1B.
For E/Z vermelhotin (6), the kinetic analysis data fit for a competitive parabolic inhibition; the statistical analysis of this fit predicted that two molecules of 6 bound to the active site (Table 2), which is possible due to the small molecular size.Figure 4A presents the mechanism of inhibition for 6, implying mutually exclusive binding of the inhibitor and the substrate, allowing for two molecules of the inhibitor to bind to the enzyme.In the doublereciprocal plot, the V MAX value did not show any change, but the K M value changed significantly (Figure 4B).The modification of the structure of E/Z vermelhotin could improve its primary site of interaction.
The classic Michaelis-Menten plots were also convenient for defining the mode of inhibition of compounds 1, 3, and 6 and have been included in Supplementary Figure S3.

SAS map
A SAS map (Medina-Franco, 2012) was created to identify scaffold hopes (compounds with low structural similarity but low activity difference) by plotting the absolute value of the PTP1B pIC 50 difference (|ΔActivity|) against ursolic acid structural similarity (Tanimoto similarity indexes) (Figure 1 and Figure 5) employing atom-type fingerprints.The resulting map was divided into four zones (ZI-ZIV), comprising the scaffold hopping region (ZI), smooth or continuous SAR (ZII), a non-descript region (ZIII), and activity cliffs (ZIV) (Naveja et al., 2018).
The analysis of the ZI region (818 pairs) revealed that NPs 3, 5, 6 (Figure 5B), 7, 8, 9, 11, 34, 82, and 99 (Supplementary Figure S15) possessed low structural similarity to the reference compound (UA) and good PTP1B inhibitory potency, therefore they are considered scaffold hops.In this group, it is not possible to make conclusions regarding structural activity relationships due to the diversity of structures.The structural simplicity of compounds 3 and 6 makes them hit compounds, likely to be optimized and used for developing potential drugs with novel scaffolds.
The SAS map allows the identification of canophyllol (1), also a pentacyclic triterpenoid, as the only true analog of UA (located in ZII), justified as the smallest ΔA obtained (0.2); the results suggest that the pentacyclic structure and oxygenated functionalities at C-3 and C-28 make similar compounds more active.This information is consistent with other UA derivatives possessing alcohol and carboxyl functionality at C-3 and C-28, respectively (Guzmán-Ávila et al., 2018;Khwaza et al., 2020).
This research also examined several triterpenoids that included 4, 10, 56, 64, 85, 86, and 89, with a lower activity level than UA.None of them had the same pentacyclic core as UA and lacked the alcohol or carboxylic functionality at C-3 and C-28, respectively, suggesting that these structural features are important for activity (Figure 6).

Molecular docking analysis of hit compounds
To predict the binding modes between compounds 1-11 and PTP1B, a docking analysis was performed using the AlphaFold enzyme (code: Q9PT91), which was first subjected to molecular dynamics simulation (MDS) (Supplementary Figure S5A,B).After 100 ns, the region between amino acids 300 and 400 of PTP1B (Supplementary Figure S6A), which is an intrinsically unstructured site, was more stable (Supplementary Figure S6B).The robustness of the structural model was supported by Ramachandran plots and quality scores, before and after the MDS.Therefore, the model shown in Figure F6B was used for docking studies.This section shows only the results of 1, 3, and 6 (Table 3; Figures 7, 8) to correlate the findings with those of the kinetic analysis.Compounds 1 and 3 bind in a pocket formed by the 3α, 6α, and 7α helices near the C-terminal region of PTP1B, as previously reported for UA (Liu et al., 2022); this outcome is in agreement with the non-competitive mechanism found in the kinetic studies.Compounds 1 and 3 interacted with different amino acids, but the nature of the interactions was predominantly hydrophobic.Compound 1 had lower energy binding.
Since compound 6 undergoes interconversion between the E (6a1, 6b1) and Z (6a2, 6b2) isomers (Supplementary Figure S4), forming an equilibrium E/Z mixture with a ratio of 6.4: 3.6 (Leyte-Lugo et al., 2012), the docking analysis was undertaken with the four more stable conformational isomers.For all conformers, 6a1, 6a2, 6b1, and 6b2, the low-energy docking events fell in two contiguous relatively separate pockets designated as I and II (Figure 8 and Supplementary Figure S7-S10), where simultaneous binding of two molecules might occur.Conformers 6b1 and 6b2 predominantly bind to pocket II, which includes active site residues supporting the parabolic competitive inhibition kinetic mechanism found for NP 6. Pocket I, where 6a1 and 6a2 might attach, contains secondary residues different from the active site (Liu et al., 2022).The lowest energy pose was in pocket II (active site), followed by the pose at site I, which may be considered an allosteric spot (Table 3; Figure 8, Supplementary Figure S7-S10).This outcome is relevant because the bidentate PTP1B inhibitors bind simultaneously to the catalytic and allosteric sites, conferring higher specificity and potency to these inhibitors (Chen et al., 2018;Akyol and Kilic, 2021).
The docking studies for compounds 2, 4, 5, and 7-11 are given in Supplementary Figure S11, S12 and Supplementary Table S3, S4.All compounds bound to an allosteric area in the enzyme; 2, 4, 5, and 7-11 presented ΔG between −7.2 and −5.6 kcal mol −1 , and their interactions were hydrophobic; metabolites 4, 7, 9, and 10 attached in the region formed by helices 3α, 6α, and 7α, the same as UA.In the case of 2, 8, and 11, they were positioned near the 3α and 9α helices and two β strands.Finally, among all hits, flavonoid 8 showed the best binding energy to the protein.

Medicinal chemistry and ADMET predictions of most active NPs
To estimate the safety and efficacy of the 11 active NPs, their physicochemical, pharmacokinetic, and pharmacodynamic profiles were analyzed to predict their potential use in drug development (Supplementary Figure S13, Supplementary Table S5-S11).All the results predicted were only taken as an indicator of passive absorption properties because NPs or NP-derived drugs have more than two violations of the Lipinski's rule and their physicochemical, pharmacokinetic, and pharmacodynamic profiles are different from those of the approved drugs (Boufridi and Quinn, 2018).
The bioavailability radar (Figure 9) predicted that compounds 1-11 possess good physicochemical properties.However, the degree of unsaturation and solubility values were out of the optimal range in almost all the molecules but 7. Regarding the medicinal chemistry rules, only 3, 6, 7, 9, and 10 fulfilled the Lipinski and Golden Triangle rules (Supplementary Table S6); compounds 2, 5, and 8 are panassay interference substances (PAINS) due to the catechol fragment motifs; in consequence, they could show false positive results in biological assays.Finally, all compounds, excluding 1, showed one or two BRENK alerts, so they possess unwanted fragments (Supplementary Table S6).
The predictions related to passive permeability in the GI tract and BBB, as well as interaction with the glycoprotein P (PGP), are presented in the BOILED-Egg (Figure 10); the visual analysis of the BOILED-Egg graphic shows that 3, 9, and 11 could passively cross the BBB, therefore they could be neurotoxic.It is possible that 5-7 and 10 may experience passive gastrointestinal absorption because they share similar physicochemical properties with oral medications, such as M. T and E fell outside the graphic, precluding any prediction regarding their absorption type.Finally, only the prenylated xanthone 7 and triterpenoid 10 could be a substrate of PGP, indicating low possibility for unsafe drug-drug interaction and problems with their excretion.
In addition, the medium lethal concentration doses (LD 50 s) were estimated using TEST v.5.1; all calculations are given in Supplementary Table S10, S11.The toxicity conclusion must be interpreted cautiously; no compound is free of toxic effects.According to these results, compounds 1, 6, and 10 were the NPs with fewer toxicological alerts.The remaining compounds presented one or two alerts.The most relevant toxicological alerts were for compounds 9 (high mutagenic and teratogenic effects), 5, and 7 (high mutagenicity), and compound 2 (high teratogenic action).Regarding the LD 50 s, compounds 1, 3, and 6 showed LD 50 s ≥ 1,000 mg/kg.
Figure 11 and Figure 12 display six toxicological properties on a graph to better understand potentially toxic molecules.On average, NPs 1 and 6 showed the lowest number of alerts (three and six, respectively) (Supplementary Table S11).Experimental Lorke toxicity evaluation of compound 6 revealed no acute toxic effect on mice; the experimental LC 50 value was estimated to be higher than 5,000 mg/kg, in agreement with the predicted results.

Conclusion
From a small and diverse NP library selected mostly on the basis of an ethnopharmacological criterion, 11 hit compounds (1-11) were identified, seven of which were scaffold hops (3, 5, 6, 7, 8, 9, and 11).NPs 1-5, 8, 10, and 11 are used in different traditional preparations for treating diabetes or its complications.Therefore, this work provides information supporting the rational use of some Mexican medicinal plants employed for treating diabetes.
The 47 active compounds detected after the preliminary screen are chemically diverse, covering most of the explored space of Mexican NPs.The Tanimoto similarity indexes and the absolute value of the PTP1B inhibitory activity differences of all pairs of compounds allowed identifying 10 hit molecules (3, 5, 6, 7, 8, 9, 11, 34, 82, and 99) with low structural similarity to the reference compound (UA) but good PTP1B inhibitory potency (scaffold hops).The structural simplicity of 3 and 6 pinpoint them as hits for developing potential drugs with novel scaffolds.
After thoroughly analyzing the medicinal chemistry and physicochemical, pharmacokinetic, and toxicological properties of various NPs, 1, 6, and 10 seem to be the most promising hit compounds as PTP1B inhibitors.According to the Lorke analysis, compound 6 is not toxic to mice, with an LC 50 value higher than 5,000 mg/kg.The three compounds warrant drug optimization analysis to develop new leads for drug development.

FIGURE 1
FIGURE 1Chemical space 3D representation: the black spots are the most active compounds tested, and gray spots represent the Mexican natural products database BIOFACQUIM.

FIGURE 3
FIGURE 3 Double-reciprocal plots of PTP1B inhibition at different concentrations of compounds 1 and 3.

FIGURE 2
FIGURE 2Compounds showing pIC 50 in the same magnitude order as that of ursolic acid (UA).

FIGURE 4
FIGURE 4 Study of kinetic inhibition for 6. (A) Equilibria for the parabolic competitive inhibition.(B) Inhibition patterns of PTP1B fitted to a global total competitive parabolic inhibition.
FIGURE 5 (A) SAS map shows the number of compound pairs in each zone.(B) Chemical structures identified as key scaffold hops (ZI).(C) Activity cliff compounds located in ZIV.

FIGURE 7
FIGURE 7Predicted interactions between ligands and PTP1B (AlphaFold code: Q9PT91).The protein is displayed in blue (cartoon), and active sites are shown as yellow spheres and compounds as sticks: red (1), green (3), and UA (orange).

FIGURE 9
FIGURE 9Visual representation of the optimal range for each property presented in axes.Trodusquemine (T), ertiprotafib (E), and metformin (M).

FIGURE 11
FIGURE 11Graphic representation of some of the estimated toxicological properties for compounds 1-11.Correlation plot of compounds according to their estimated LD 50 and toxicity class.Marker size indicates the probability for a compound to be a PGP inhibitor (the larger, the more probable).Marker color refers to the toxicity score (chemical alerts found in the structure; values > 1 indicate unfavorable substructures).Marker shape indicates high (circle), low (square), or negligible (triangle) mutagenicity.Marker background color correlates with the probability of being a PAINS.

FIGURE 12
FIGURE 12Graphical representation of some of the estimated toxicological properties for compounds 1-11.Correlation plot of compounds according to their estimated LD 50 and toxicity class.Marker size indicates the Caco-2 cell permeability related to in vivo absorption.Marker color refers to low or nontumorigenic effect.The marker background color indicates teratogenicity (reproductive effect).

TABLE 1
In vitro inhibitory activity against PTP1B of most active compounds of the natural product library.